Competition of hydrophobic and Coulombic interactions between nano-sized solutes 
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The solvation of charged, nanometer-sized spherical solutes in water, and the effective, solvent- 
induced force between two such solutes are investigated by constant temperature and pressure 
Molecular Dynamics simulations of model solutes carrying various charge patterns. The results for 
neutral solutes agree well with earlier findings, and with predictions of simple macroscopic consid- 
erations: substantial hydrophobic attraction may be traced back to strong depletion ("drying") of 
the solvent between the solutes. This hydrophobic attraction is strongly reduced when the solutes 
are uniformly charged, and the total force becomes repulsive at sufficiently high charge; there is a 
significant asymmetry between anionic and cationic solute pairs, the latter experiencing a lesser hy- 
drophobic attraction. The situation becomes more complex when the solutes carry discrete (rather 
than uniform) charge patterns. Due to antagonistic effects of the resulting hydrophilic and hy- 
drophobic "patches" on the solvent molecules, water is once more significantly depleted around the 
solutes, and the effective interaction reverts to being mainly attractive, despite the direct electro- 
static repulsion between solutes. Examination of a highly coarse-grained configurational probability 
density shows that the relative orientation of the two solutes is very different in explicit solvent, 
compared to the prediction of the crude implicit solvent representation. The present study strongly 
suggests that a realistic modeling of the charge distribution on the surface of globular proteins, as 
well as the molecular treatment of water are essential prerequisites for any reliable study of protein 
aggregation. 
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I. INTRODUCTION 

It is a well-known fact of physical chemistry that solvo- 
phobic solutes of similar sizes and shapes tend to attract 
each other in an incompatible solvent. Classic examples 
are the effective attraction between monomers of poly- 
mer coils in poor solvent, which leads to collapse below 
the 0-temperature^ or the attraction between hydropho- 
bic surfaces in water £ The effective attraction ultimately 
leads to phase separation of the solvent and solute as 
the concentration of the latter increases. The solvent- 
averaged effective interaction (or potential of mean force) 
is related to the variation of free energy upon bringing 
the two solutes from infinite to a finite separation in the 
solvent. The change in free energy has an entropic com- 
ponent, associated with the reorganization of the solvent 
molecules around the two solutes, and an energetic con- 
tribution which accounts for the deficit in attractive in- 
teractions between solvent molecules close to the solutes. 
There is some analogy between solvophobic attraction 
and the well-known depletion interaction between col- 
loidal particles induced by a depletant like non-adsorbing 
polymers^ In the latter case the depletion attraction can 
be essentially understood in terms of excluded volume, 
and is hence of entropic origin, while hydrophobic in- 
teractions have a large energetic contribution, associated 
with the formation or break up of hydrogen bonds. 

It has been recognized that the size of the solute plays 
an important role in understanding its solvation energy, 
and effective solute-solute attraction^ For solutes of a 
characteristic size larger than a few molecular diame- 
ters (typically larger than lnm), a mechanism first en- 
visioned by Stilinger— is that of solvent dewetting ( "dry- 
ing"), i.e. the solvent molecules tend to move away from 



the surface of a large solute, and form a liquid-gas like 
interface parallel to the solute interface The overlap 
of the drying zones associated with two large solutes as 
their surfaces come together may then give rise to an ef- 
fective attraction^ very much like the depletion mecha- 
nism between-colloidal particles. The mechanism holds, 
a priori, for any solvent, provided that it is in a ther- 
modynamic state close to liquid-vapor coexistence^ The 
drying mechanism and resulting attraction between two 
plate-like solutes was confirmed by Molecular Dynamics 
(MD) simulation of Wallquist and Berne. 8 Similar simu- 
lations were carried out for two large spherical solutes in 
a Lennard- Jones solvent, and attractive solvation force 
profiles were determinedAifliii 

However most nano-scale biomolecular solutes, like 
proteins, carry electric charges, which make them, at 
least partly, hydrophilic. The main objective of the 
present paper is to investigate the influence of solute- 
solute and solute-solvent electrostatic interactions on the 
effective, solvent-induced potential of mean force between 
two solutes. In order to make contact with earlier work on 
neutral hard-sphere solutes, we restrict the present inves- 
tigation to spherical solutes of identical radii R < lnm, 
but carrying various surface charge patterns. The two 
main questions which will be addressed are: a) how does 
the competition between hydrophobicity and electrostat- 
ics affect the total effective force between anionic and 
cationic solutes; b) is the total force sensitive to details 
of the charge patterns carried by the solutes? In par- 
ticular are there significant differences between results 
obtained with continuous and discrete charge patterns of 
the solutes? Such differences were recently highlighted 
by a calculation of the second virial coefficient in an im- 
plicit solvent model of globular proteins, which does not, 
of course, allow for hydrophobic attraction^ 
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FIG. 1: Sketch of the SPC/E water model. The 
vector Co embodies the orientation of the molecule. 



We have attempted to answer these questions by a se- 
ries of constant pressure MD simulations of two spherical 
solutes of varying radii (up to R = 1.3nm) immersed in 
water modeled by the SPC/E intermolecular potentials- 
taken under normal conditions, i.e. close to liquid-vapor 
coexistence. The paper is structured as follows. The 
models and simulation procedures are detailed in Sec. 
ITT1 The solvation of single charged solutes is examined 
in Sec. IHII The effective interaction between two solutes 
as a function of the mutual distance is estimated from a 
simple macroscopic theory in Sec. IIVI where the method 
for extracting the mean effective force from simulations is 
also defined. The results from MD simulations for several 
charge patterns are presented in Sec. ]V\ while concluding 
remarks are made in Sec. IVII 

Part of the present results were briefly reported 
elsewhere^ 



II. MODELS AND METHODOLOGY 

The MD simulations were carried out on periodic 
samples containing N w water molecules and one or 
two solutes. The SPC/E model of a water molecule 13 
is sketched in Fig. 1. Two water molecules interact 
via a Lennard-Jones potential between the oxygen (O) 
sites, and the bare Coulomb potentials between the 9 
pairs of sites. The Lennard-Jones parameters are e = 
0.6502kJmor 1 and a = 3.169A. The molecules are as- 
sumed to be rigid (with OH bond lengths and HOH bond 
angle specified in Fig. 1) and nonpolarizable. The solutes 
are smooth spheres of bare radius Rq, which interact with 
the O-site of the water molecules by the purely repulsive 
potential 

V (r) =0(r-i? o )- 12 , (1) 

where r is the distance from the solute center to the O- 
site, and the energy scale 4> is chosen such that the O- 
atom experiences a repulsive energy k^T at a distance 
r — Rq =1A from the solute surface. With this conven- 
tion, the effective radius of the solutes may be defined as 
R = Rq + lA. The purely repulsive interaction Q is cho- 
sen to mimic a strongly hydrophobic interaction between 
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TABLE I: Characteristics of the different models 
So±, To±, Co± used. N c is the number of charges carried 
by the solute, while Q is the net charge. Rq is the bare so- 
lute radius. By R c we denote the distance of the charges 
to the center of the solute, and R cc quantifies the nearest 
neighbor distance between the charges. In the charged S- 
models the charge is located in the center of the solutes, 
while in the T and C models the charges are distributed 
tetrahedrally and cubically on the sphere surface. 



the neutral solute and the water molecules. The model 
involving spherical solutes with no charged site will be 
referred to as So- 

Since the main objective of our work is to investigate 
the difference in the effective, solvent induced interac- 
tion between the cases of neutral and charged solutes, we 
have considered several models for the latter (cf. Table 
1). In the simplest model, a total charge Q = qe (where 
e is the proton charge) is assumed to be uniformly dis- 
tributed over the solute surface. According to Gauss' 
theorem, this is equivalent to placing a single charged 
site Q at the center of the spherical solute. We consider 
both anionic (q < 0) and cationic (q > 0) solutes and 
the corresponding models will be referred to as S_ and 
S+ . To ensure overall electroneutrality of the system, the 
total charge carried by the solutes must be compensated 
either by a uniform background of total opposite charge 
permeating the system, or by explicitly including coun- 
terions. For the latter we choose Cl~ (cationic solutes) 
and Na + (anionic solutes) ions. Their mutual interac- 
tions and coupling to water molecules involve the stan- 
dard Coulomb-interactions, and a Lennard-Jones part, 
with e and a parameters taken from Spohr^ The short 
range interaction with the solutes is again described by 

CD- 

In order to investigate the sensitivity of the effective 
forces to details of the solute charge patterns, we have 
also considered models with discrete charge distributions 
involving N c point charges placed on the surface of the 
solutes. In the tetrahedron model, N c — 4 charges are 
tetrahedrally arranged at a distance Rq from the cen- 
ter of the solutes, as illustrated in Fig. |2{a) ; we consider 
both the cases where all 4 charges are of the same sign 
(q = ±4, referred to as models T + and T_), and the neu- 
tral situation, where two charges are positive and two 
are negative (model To). We have also considered cu- 
bic charge distributions, as sketched in Fig. Etb) , where 
N c = 8. In models C+ and C_ all 8 charges are positive 
{q = 8) or negative (q — —8), whereas in model Co, four 



3 




FIG. 2: Solute with (a) tetrahedral (T mod- 
els) and (b) cubic (C models) charge distribu- 
tion. The charges are connected in this fig- 
ure for a better visualization of the structure. 



vertices carry a charge +e, while the other four carry op- 
posite charges in an alternating arrangement such that 
the three nearest neighbors of a negative charge are posi- 
tive and vice versa. A similar model of globular proteins 
was considered by Allahyarov et al.^ but in an implicit 
(continuous) solvent representation. 

In order to avoid very close approaches of water H-sites 
and the solute surface charges, the charged sites at the 
vertices of the tetrahedron or cube, situated at a distance 
Rq from the solute center, are not simply point charges, 
but are modeled by Cl~ or Na + ions. The corresponding 
LJ potentials prevent these sites and the water H atoms 
to come too close, and hence unreasonably large electro- 



static forces, which could lead to electrostatic "sticking" 
of the water molecules to the solute surface. The total 
interaction energy between one solute and the N w water 
molecules in a periodically repeated, cubic simulation cell 
is: 

v sol = Ew+EE^? 1 ) 

i—l i—l a—1 

N w N c 3 

+ EEE^^w(rf), (2) 

8=1 a=l 0-1 

where r 8 - is the distance from the center of the solute to 
the O-atom of the ith water molecule, and r" is the dis- 
tance from site a on the solute to site [3 of the ith water 
molecule ((3=1 for the oxygen site). The first term on 
the rhs of Eq. J2J) corresponds to the short ranged repul- 
sion QJ; the second term is the sum of Lennard- Jones 
interactions between all N c sites of the solute and the O- 
sites (/? = 1) of the N w water molecules, which depend on 
the corresponding site-site distance rf 1 ; finally the last 
term accounts for the Coulombic interactions between all 
N c solute sites and all 3 water sites; (j>Ew(r) is the electro- 
static interaction between two elementary charges, prop- 
erly summed over an infinite array of periodic images, us- 
ing the smooth-particle-mesh Ewald (SPME) methodic 
(see the Appendix A for details). An expression similar 
to Eq. (0) holds for the total interaction between a solute 
and its CI - or Na + counterions. 

The MD simulations were carried out with the 
DLPOLY2 17 package, using the Verlet leapfrog 
algorithm ji£ with a timestep of 2fs. Simulations 
were carried out at constant pressure (P = latm) and 
constant temperature (T = 300K), using appropriate 
barostats and thermostats (see Appendix A). We 
emphasize the importance of using constant pressure 
simulations of charged, aqueous systems: electrostriction 
and drying mechanisms modify the density of water in a 
finite, closed system in a significant way. In the present 
NPT ensemble simulations, the overall density of water 
varied from p = 0.033A to 0.035A" ; in estimating 
the average water density, the effective volume AttR 3 /3 
of the solutes must be subtracted. The choice of the box 
length L of the periodically repeated simulation cell and 
the arrangement of solutes inside the cell are discussed 
in the Appendix A. The cell contained up to A^ w = 3000 
water molecules. 



III. SOLVATION OF A SINGLE CHARGED 
SOLUTE 

Before embarking on the main subject of this paper, 
namely the water-induced effective interaction between 
two protein-like solutes, we first consider the solvation of 
a single, neutral or charged spherical solute, a problem 
which has been abundantly addressed in the literature, 
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FIG. 3: Density profiles of water oxygen and hydro- 
gen atoms (inset) around one isolated neutral solute 
(S ) plotted for different radii R/A = 2,3,4,5,6,11. 




FIG. 4: Orientation parameter P(r) defined in Eq. 
of the water particles around a solute with radius R = 
11 A carrying no charge (So, solid line), and charge q = 
±5 (S±, dashed lines), and q = ±10 (S±, dotted lines). 



since the pioneering work of BorniS on solvation free en- 
ergies of ions, and of Reiss et al. on the scaled-particle 
theory of cavity formation and hard sphere solutesiSS 



A. Water structure around a solute 

Consider first the structure of water around an isolated 
neutral solute (model So). Fig. shows MD results for 
the oxygen and hydrogen radial distribution functions 
(RDF), or density profiles, around a solute for various 
solute radii R. The height of the main peak in the RDF of 
both O and H sites of the water molecules is seen to first 
grow with increasing R, and this may be rationalized in 
terms of enhanced packing of the molecules at the surface 
of the solute. But for R > 5A, the peak height is seen to 
decrease monotonically, due to the unbalanced attraction 
experienced by the water molecules near the surface from 
bulk water. Very similar predictions of the depletion of 
water around large spherical solutes (or cavities) have 
been reported earlier in the literatureA^ The hydrogen 
and oxygen peaks are located at nearly the same distance 
from the solute for any given radius R, with a tendency 
of the hydrogen peak to be slightly further out. This 
seems to indicate that there is no strong orientation of 
the water molecules in the first solvation shell towards or 
away from the solute surface. 

This observation may be quantified by considering the 
following orientational order parameter: 



P(r) = 



MM 



(3) 



where uj is the water molecule orientation vector (cf. 
Fig. 1), and the configurational average is taken for a 
fixed distance r from the solute center to the O-site of 
the water molecules. MD results for a neutral solute 
of radius R = llA are plotted in Fig. The curve 




FIG. 5: Density profiles of water oxygen and hydro- 
gen atoms (inset) around a (a) positively (S+) and 
(b) negatively (S_) charged solute with radius R = 
11 A and central charge q — ±5 (dashed line), and 
q = ±10 (dotted line). We also plot the result 
for the Sq model (solid line) with the same radius. 
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P(r) takes slightly positive values (P ~ 0.05) for water 
in the first solvation shell (r ~ 12. 5 A), indicating a weak 
tendency of the hydrogen atoms to point away from the 
solute. 

The effect of charging the solute is illustrated in 
Fig. ETa) and (b), where the RDFs are plotted for a fixed 
radius R — llA, and charges q — 0, ±5, and ±10, within 
the So and S± models (neutral, or uniformly charged so- 
lutes). For positive charges, Fig. (Sfa), the height of the 
first peak in both oxygen and hydrogen RDFs is seen 
to shift to shorter distances, and to increase as q in- 
creases, signaling an effective attraction of the dipolar 
solvent molecules due to the radial electric field. The 
trend is similar for negative charges, as regards the oxy- 
gen RDF. However the initial first peak (at r ~ 12A for 
q = 0) in the hydrogen RDF is seen to split into a pre- 
peak around r ~ 10 A and a broad feature close to the 
initial peak when q = —5. This points to a reorienta- 
tion of the water molecules in the first hydration shell, 
with the positive hydrogen atom preferring being closer 
to the surface of the anionic solute. Further decrease of 
the negative charge (q = —10) consolidates this struc- 
ture, with two hydrogen peaks growing in amplitude (cf. 
Fig. OJb)). Interestingly, while the hydrogen RDFs are 
very sensitive to the sign of the solute charge (S+ ver- 
sus S_), the amplitudes of the first peaks of the oxygen 
RDFs are nearly independent of this sign, but the peak 
around the positive solute appears to be broader, sig- 
naling a larger water coordination number in the first 
solvation shell. The corresponding orientational order 
parameter P(r) is plotted in Fig. 0] The absolute value 
of P{r) has maxima at contact r ~ 11 A and approx- 
imately one water diameter further away (r ~ 13. 5A), 
and increases with absolute charge, irrespective of the 
sign of the charge carried by the solute. Closer inspec- 
tion of the curves in Fig. ^reveals, however, a significant 
asymmetry, if not in the overall shape of P(r), at least 
in the amplitudes, which are typically 20% larger for the 
negative solute. Anion/cation hydration asymmetry had 
already been reported for microscopic ions in aqueous 
solutioniSi 

We now turn to the structure of water around a solute 
with an inhomogeneous charge distribution, restricting 
the discussion to the tetrahedral T + and T_ models. In 
order to characterize the anisotropy of the problem, it is 
desirable to distinguish between water molecules close 
to the four surface charges, and the remaining water 
surrounding the solute. We achieve this by averaging 
over water molecules whose centers fall either inside or 
outside well-defined cones whose axes coincide with the 
radii joining the solute center and the surface charges and 
whose vertices coincide with the solute center. A sketch 
of the two-dimensional projection of one of the sides of 
the tetrahedron and its associated cones is shown in the 
inset to Fig. EJa) . Averages are taken over cones of open- 
ing angle O = 30° , high enough to accommodate the first 
two solvation shells around a surface charge. The results 
for the oxygen density profiles for water molecules in- 
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FIG. 6: Oxygen density profiles around the tetrahedral 
(a) positive T + and (b) negative T_ solutes. The curves 
are for a full angular average (long dashed lines), cone 
averages around the charges (circles), and averages of 
water excluded by these cones (squares), as explained in 
section All Al We also plot the density profile around a 
neutral solute So (solid line) for comparison. The in- 
set in (a) sketches a two dimensional projection of a 
tetrahedral solute. The cone averages are performed 
over the water molecules in the grey cones (opening an- 
gle = 30°) containing a surface charge (dark circle). 



side or outside the cones are shown in Figs. HJa) and 
(b) for the T+ and T_ models, respectively. The water 
molecules inside the cones exhibit a typical solvation shell 
structure with a large first peak, and a much lower sec- 
ond peak, which is hardly visible in the T + case. Outside 
the four cones, water appears to be highly depleted com- 
pared to its distribution around a neutral So solute, up 
to a radial distance r ~ 13A. Thus the 25% of the solute 
surface area inside the cones act as hydrophilic "patches" 
while the remaining 75% are hydrophobic. Interestingly, 
if an angular average is taken over the total solute area, 
the mean density of water close to the surface (< 13A) of 
a T + or T_ solute is significantly smaller than the corre- 
sponding density around a neutral So solute ! Integration 
of the water density profile up to r = 13A yields coordi- 
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FIG. 7: Orientation profiles P(r) around the tetra- 
hedral positive T + (circles) and negative T_ (no 
symbols) solutes. Cone averages around charges 
(solid lines) are compared to averages of water 
molecules excluded by those cones (long dashed lines). 




FIG. 8: Solvation free energy A/zo of a neutral 
spherical solute So in water. The inset shows A/zo 
divided by the sphere surface showing an asymp- 
totic approach for large R to the liquid-vapor sur- 
face tension of water (dashed line), 7 « 72m J/m 2 . 



nation numbers (numbers of water molecules) of 133 for 
S , 92 for T+, and 95 for T_, showing a 30% depletion of 
water around the T-solutes compared to So- The orien- 
tational order parameter (3) for the T + and T_ models 
are plotted versus r in Fig [7| The orientational order 
of water molecules inside the cones is seen to be similar 
to that around homogeneously charged solutes S + or S_ 
(cf . Fig. . In the depleted volumes outside the cones 
the water molecules shows little orientational order; if 
anything they tend to orient in the direction opposite to 
the mean orientation inside the cones. 

Qualitatively similar observations hold for the distribu- 
tion of water molecules around solutes with cubic charge 
distribution (models C + or C_), but obviously the vol- 
umes depleted of water are now smaller, since the "hy- 
drophobic patches" have shrunk now to only half of the 
solute surface area. The MD simulations for the so- 
lutes with non- vanishing net charge were carried out with 
explicit counterions. Test runs where these counteri- 
ons were replaced by a uniform neutralizing background 
showed no differences within the statistical uncertainties. 



B. Solvation free energy 

The solvation free energy is equal to the reversible work 
required for transferring a solute from vacuum into a sol- 
vent. For neutral solutes in water, the solvation free en- 
ergy is generally positivefS* 2 ^ and for atomic-size solutes, 
it stems mainly from the entropy cost of the restructur- 
ing water molecules around the solute. For larger spher- 
ical solutes, a cross over in the variation of the solva- 
tion free energy with radius occurs typically around 1 
nmpi The classic Born model 19 provides the simplest ap- 
proach to the solvation free energy of charged spherical 



solutes. The solvent is treated as a dielectric continuum 
of permittivity e and the hydration free energy increases 
quadratically with solute charge and is proportional to 
the inverse of the Born radius Rb, according to 

-iSs* 1 - 1 **- (4) 

Note that the solvation free energy is a difference in chem- 
ical potential of the solute as it is moved from vacuum 
into the solvent. Hydration free energies from the Born 
model agree well with experimental values, once the un- 
known parameter Rb is defined. The Born radius for an 
ion can deviate substantially from its Pauling radius (a 
measure of the size of an ion)i2i 

We have obtained solvation free energies for our model 
solutes by thermodynamic integration, using the general 
formula 

where the coupling parameter A gradually "switches on" 
the interaction (2) between the solute and the solvent 
from an initial state A — Xq (say a neutral point solute) 
to a final state A = Ai corresponding to the complete so- 
lute/solvent system. The brackets < .. >a denote a sta- 
tistical average over all solute-solvent configurations for 
a solute-solvent coupling characterized by V so i(A). The 
index sim in A/i s ; m indicates that the estimate of the sol- 
vation free energy is based on MD simulations of a finite 
sample; finite size corrections will be added as explained 
later. 

In practice, we proceded in two steps. In a first stage, 
we computed the solvation free energy A/xo of a neutral 
spherical solute (model So), as a function of its radius 
R. The second step is to charge up the initially neutral 
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FIG. 9: Excess solvation free energy Afj,± of charging a 
spherical solute of radius R — llA in water homoge- 
neously to a charge q (solid line) and — q (long dashed 
line). The symbols denote the solvation free energies of 
charging the To± and Co± models in water. The symbols 
at q = 4 are for the tetrahedron, at q ~ 8 for the cube. 
The overall charge of the T ± and Co± solutes is zero 
(circles), positive (diamonds), and negative (squares). 



solute to the final charge pattern. In step one, the cou- 
pling parameter A in Eq. JSJ) is simply the radius itself, 
and A/x is consequently the work required to blow up 
the solute against the normal force exerted by the solvent, 
integrated over the particle surface; in this case the force 
is just the radial derivative of the first term on the rhs 
of eq (J2J). This is implemented, in practice, by starting 
from Ao = R = 0, and increasing the radius by steps of 
AA = AR = lA, up to Ai = 14A (the largest neutral so- 
lute considered in the present work). The averaged radial 
force, as obtained from the MD simulations for various 
R, is interpolated with a cubic spline and integrated to 
yield A/iq. The resulting solvation free energy is plotted 
in Fig. |SJ and is seen to increase monotonically with R. 
The solvation free energy per unit area, A^o/4ttR 2 , is 
plotted in the inset to Fig. |SJ and is seen to approach 
asymptotically a constant value for radii R > 10A; the 
latter is close to the liquid-vapor surface tension of wa- 
ter, 7 = 72m J/m 2 . This behavior is close to that re- 
ported by Lum et al& for hard sphere cavities in water. 
The "softer" solute-solvent pair potential (1) used in the 
present work does not modify the solvation process sig- 
nificantly, compared to the case of hard sphere solutes. 

In the second stage, to go from the neutral to the 
charged solute, the charges of the N c sites on the solute 
are gradually turned on, i.e. (? Q (A) = Xq a (1 < a < N c ), 
where A is varied from to 1. The quantity to be aver- 
aged in Eq. (JSJ is now the total electrostatic energy of the 
solute in the field of the water molecules and their peri- 
odic images; the statistical average is to be taken over all 
Boltzmann-weighted water configurations when the elec- 
trostatic solute/solvent coupling is multiplied by A. Since 



for any A > the system carries a net charge, a compen- 
sating uniform background charge must be included in 
evaluating the Coulombic part of Eq. @ by Ewald sum- 
mation. If the self interaction energy of the solute with its 
own images and the neutralizing background is properly 
included, the resulting free energies are virtually indepen- 
dent of the size L of the simulation boxjS 2 ^ see the Ap- 
pendix B. Results for Afj,± corresponding to a uniformly 
charged solute with radius R = 1 1 A are plotted in Fig. 
as a function of q, for anionic and cationic solutes. Note 
that this is the excess free energy for charging an initially 
neutral solute of R = 11 A, previously inserted into the 
solvent, to a charge q. In order to obtain the total solva- 
tion free energy, the contribution A/io corresponding to 
the insertion of the neutral solute in water (« 0.6MJ/mol 
for R = 11 A) must be added to Afx±. The curves are 
essentially quadratic on the scale shown, in agreement 
with Born theory. Least squares fits of the data to the 
Born formula (4) yield R^ = 9. 8 A and R^ = 8. 6 A , both 
smaller than the effective solute radius R = 11 A. The sol- 
vation free energy of cationic solutes is slightly positive 
for < q < 1 . Such a behavior has already been reported 
for small cationic solutes (e.g Na + )22i2i and may be un- 
derstood from the competition between the free energy 
cost of the rearrangement of water around the solute, and 
the electrostatic energy gain of the dipolar solvent in the 
electric field of the solute. The latter contribution ap- 
pears to dominate already for small \q\ in the case of an- 
ionic solutes, for which A/i_ is always negative. Over the 
whole range of absolute charge A/x_ is systematically 
lower than A/i + , pointing to a preferential solvation of 
anionic solutes, again in agreement with earlier findings 
for other charged solute mo dels S 2 ^ 2 ^ In Sec. IIII AI we 
learned that the restructuring of water is stronger around 
an S_ solute than around its S+ counterpart, so that one 
would expect a higher cost in entropy. Apparently the 
closer approach of the hydrogen atoms to the negative 
solute decreases the electrostatic contribution to the free 
energy more than the positive entropic cost, resulting in 
an overall lower solvation free energy. 

Solvation free energies for solutes carrying discrete 
tetrahedral or cubic charge distributions, corresponding 
to models T+, T_, and To, and C+, C_, and Co are also 
shown in Fig. |5| The contribution to the solvation free 
energy from the steric LJ-part of the surface charge inter- 
action with the water molecules is insignificant and was 
neglected in the calculation of A/j,±. As in the case of the 
uniformly charged solutes, the negative solutes are prefer- 
entially solvated compared to their positive counterparts. 
The solvation energies of the overall neutral solutes To 
and Co lie well above those of the charged solutes (T± or 
C±). In particular \Afi\ of the overall neutral solute with 
a cubic charge pattern Co is roughly three times smaller 
than the corresponding |A/x±| of the globally charged so- 
lutes C+ and C_. This may be a consequence of the 
considerable reorganization of water around the solutes 
with surface charges of alternating sign, resulting in a 
substantial cost in free energy. Finally, the solvation free 
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energies of solutes with discrete charge patterns (T± or 
C±) are seen to lie 10-20% below the solvation free en- 
ergies of their uniformly charged counterparts (S± with 
q = ±4 and ±8). 



IV. EFFECTIVE INTERACTION BETWEEN 
TWO SOLUTES 

We now turn to the main objective of this pa- 
per, namely the determination of the effective, solvent- 
mediated interaction between two nanometer-sized neu- 
tral or charged solutes in water. In subsection IV A and 
IV B we derive this interaction from simple macroscopic 
considerations, while the MD methodology is presented 
in IV C. 



A. Phenomenological theory for charged plates 

A simple macroscopic argument, similar to Kelvin's 
theory of capillary condensation predicts that water near 
liquid/vapor coexistence will undergo "drying" when 
confined between two hydrophilic plates, below a critical 
distance D c separating these plates— We extended the 
argument to the case of charged plates^ showing that 
D c is strongly reduced by the electrostatic energy asso- 
ciated with the surface charge carried by the plates. The 
macroscopic argument is further refined hereafter. Con- 
sider two parallel plate-like solutes of area A\, separated 
by a distance D and carrying opposite surface charges 
±er, immersed in a polar solvent of dielectric permittiv- 
ity e. Neglecting edge effects (an approximation valid as 

1/2 

long as D C A x ' , the electric field between the plates 
is Eq/c with E = <r/e n . We require the difference in the 
grand potential between the situations where the liquid 
solvent (I) or its vapor (g) fill the volume A X D between 
the two plates: 

1 E 2 

0* = -P a A x D + -eo-MiD + 2-y wa Ai+~/ la A 2 , (6) 
z e Q 

where P a is the pressure of phase a = l,g and j wa the 
surface tension between phase a and the plate ( "wall" ) . 
A2 is the area of the liquid-vapor interface limited by 
the edges of the two opposite plates, which is created 
when the volume between the plates is filled by vapor. 
jl a is the liquid-vapor surface tension when a = g, and 
vanishes of course when a — I. The last term in Eq. © 
may be neglected for infinitely large platesJ^ Consider a 
state close to phase coexistence at temperature T, and 
let S/j, = ji— /i sa t be the positive deviation of the chemical 
potential from its saturation value. Expanding the P a to 
linear order around their common value at saturation, 
one arrives at the following expression for the difference 



in grand potentials: 

Q g -Ql = ( Pl -p g )5 f iA 1 D + ^-E 2 {---)A 1 D 

2 e g ei 

+ 2(j wg - -f w i)Ax ~ ji g A 2 . (7) 

In water e; = e ^> e g ~ 1, p g <C pi (except near critical 
conditions) and j w i — j wg — 7; g = 7 for a purely hy- 
drophobic surface. Moreover A2 = UD, where U is the 
circumference of one plate. Hence: 

Ail = (piSpi + ^£ 2 ) A X D - 2 1 A 1 + jUD. (8) 

AS1 = fig — Cli is the reversible work bringing the two 
plates from infinite separation (when the volume between 
them is filled by liquid) to a distance D, at which "dry- 
ing" has already occurred. At contact, Afl(D = 0) = 
— 2jAi, and Afi increases linearly with D. The range of 
the purely attractive potential is defined by the distance 
D c at which Af2 = 0; for D > D c , the liquid is the pre- 
ferred phase between the plates; for D < D c , "drying" 
occurs. From Eq. (JSJ we obtain 

D 2 -2 (9) 

c " ptSfi + ^E 2 + r/U/Ax ■ 1 ' 

Near the liquid-vapor transition of water Sfi <C fcsT, and 
may be neglected compared to the surface tension term 
in the denominator. Consider circular plates of radius R; 
then U/Ai = 2R, and if they are uncharged (Eq = 0), 
D c = R, which is in agreemeent with previous calcu- 
lations of the mean force between plate-like solutes in 
water— and in a LJ-fluid. 26 In the case of high surface 
charges (a < e/nm 2 ), Eq can be as large as 10 10 V/m, 
and the corresponding electrostatic term in the denomi- 
nator becomes comparable to the surface term for solute 
sizes of a few nm. This leads to a strong reduction of 
D c compared to the case of neutral solutes. This re- 
duction of D c hints at a considerable weakening of the 
hydrophobic interaction between two solutes when the 
latter are charged. This trend will be confirmed by the 
MD results in Sec. V. Note however that the simple 
macroscopic model ignores molecular details, and that 
its prediction does not, a priori, apply to solutes carry- 
ing discrete charge patterns (i.e. "hydrophilic patches") 
for which a more microscopic description is required. 



B. Phenomenological theory for spherical solutes 

The previous model can be extended to the case of 
neutral or charged spherical solutes with radius Rq as fol- 
lows (see Fig. llOf) . When "drying" occurs, the simulation 
data (cf. Fig. IT3I (b) or (c)) suggest that a cylindrically 
symmetric domain bounded by the two spherical solute 
surfaces (Si) and the curved liquid-vapor meniscus (£2) 
is filled with vapor. The surface S2 is assumed to touch 
the solute spheres tangentially (contact angle 71") and to 



9 





\R' 








\v 






.1 . - 









7 



FIG. 10: Sketch of two spherical solutes of radius 
Rq at a surface-to-surface distance s. The white 
volume V between the solutes approximates the re- 
gion depleted of water. Si and S2 are the sur- 
rounding solute- vapor and liquid- vapor surface areas, 
resp. R' is the radius of the curved surface 15*2. 




FIG. 11: Free energy Af2 of dried states between neu- 
tral spherical solutes with radius Ro = 10 according to 
Eq. Qllfl. The volume of the empty state is described 
by the parameter x for a given geometry, see Fig. 1101 
The most bottom curve is for a surface-to-surface dis- 
tance s = 0, then s is incremented by 0.5A steps. 



have a radius of curvature R' . For a given surface-to- 
surface distance s of the two solutes, the volume V of 
the "dry" domain and the areas Si and S2 are conve- 
niently expressed in terms of the single parameter X ■ clS 
depicted in Fig.^| If x — 0, the vapor domain shrinks 
to zero, i.e. the space between the solutes is filled with 
liquid, while for x = Rq the vapor occupies a cylindrical 
volume V — ttRq[(2R + s) — 4R /3]. For intermediate 
values of x, the areas Si and 62 are given by 



Si(x) = AttRqx 



S 2 (x) 

r! 
ti 



AttR' 



h arcsin 



a/2 



R' 



(x + s/2) 



Ro(x + s/2)/(R -x) 
Ro + R 



Ro 



■V2R, 



qX — X* 



(10) 



The difference in grand potential between the "dry" 
and filled states is then given (in absence of the electric 
charges) by the following generalization of Eq. JHJ : 



An = pi5fiV(x) - j s Sx(x) + r )S 2 {x), 



(11) 



where 7 S is the surface tension of the solute-liquid in- 
terface and 7 the liquid-vapor surface tension. The first 
term in Eq. i|ll|) is the bulk free energy for creating a cav- 
ity of volume V(x) in water, and favors the filled state. 
Again, near liquid-vapor equilibrium Sfi <gC ksT, and the 
volume term may be safely neglected. The second and 
third terms in Eq. I|llf> are the surface free energies for de- 
creasing the solute-liquid and increasing the liquid-vapor 
interface, respectively. For simplicity we first assume 
7 = 7 S 7 ; effects arising from 7 / 7 S will be discussed 



later. In the following we use the liquid-vapor surface 
tension of water 7 = 0.174/cbTA . AQ(x) is plotted 
versus the geometric control parameter x in Fig. llll for a 
solute of radius Rq = lnm, and various surface-to-surface 
distances s. At contact (s = 0), Afl(x) exhibits a single 
negative minimum at the non-zero value x = x m ; n (s = 0) 
signaling that the "dry" state with volume V(x m i n ) is 
stable. As s increases, the global minimum is raised to 
less negative energy values and occurs at smaller values 
of x, while a second local minimum appears at x = 0, 
corresponding to a metastable, filled state. At the crit- 
ical value s c ~ 3.2A, the global minimum jumps from 
the non zero value of 1 to 1 = x m i n (s c ) = 0. For 
s > s c , the filled state is favored. The effective interac- 
tion w(s) between two solutes is equal to the free energy 
difference in bringing them from infinite separation to a 
surface-to-surface distance s, and is given by Eq. Ullfl . 
i.e. w(s) — AQ(s; x m i n (s)) . The energy at the global 
minimum is negative for all s < s c (x m i n > 0), so that 
the interaction is always attractive. s c is the range of the 
interaction; at s = s c , w(s) = Afl = 0, and j s Si = 7S2. 

The effective potentials and forces between two solutes 
of different radii R = 5, 8, 10, 12A from Eq. ITTJ are 
plotted in Fig. A detailed numerical investigation 
shows that, assuming 7 S = 7, the range of the potential 
(and of the resulting effective force) is s c ~ 0.32i?o, scal- 
ing linearly with solute radius, independently of 7. The 
contact value of the potential is w(0) ~ — 1.197^, scal- 
ing with the solute surface area. The contact value of the 
force is F(0) ~ — 4.287-Ro- The force increases roughly 
linearly with a slope independent of Rq (cf. Fig. 12). The 
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FIG. 12: Results from the phenomenological theory in 
section IIVBI for the mean force, F* — 0FK, between 
two neutral (q = 0) spheres in SPC/E water. Results 
are plotted for solute radii Rq = 5A , Rq = 8A , 
i?o = 10 A , and i?o = 12A . The inset shows the in- 
tegrated force (potential of mean force). The depth 
and range of the force and potential increases with R. 



results are accurately represented by 

-cliRq + a 2 s for s < s c = a^Ro; 



F(s) = 7 







otherwise. 



and 



, . I — aoR% + 0,-lRqs — for s < s c ; 

w(s) = 7 < z 

otherwise. 



(12) 



(13) 



with do = 1.19, ai — 4.28, 03 = 0.32, and s c — C13R0. 
If the force is assumed to be linear in s, a 2 is fixed by 
the other constants according to a 2 = 2(0103 — aa)/a\ ~ 
3.51. The linear scaling of the contact value of the force 
with solute radius R may be rationalized by a simple 
consideration of the potential of mean force for plates 
(R = 00) at contact, w(0) — —2j and an application 
of the Derjaguin approximation, valid for weakly curved 
substrates (i.e. large R)^ this leads to the estimate 
F(s — 0) = — 2irR"f, which indeed predicts linear scaling. 

Near a "realistic" solute, water will have a surface ten- 
sion 7 S ^ 7, due to Van-der-Waals and electrostatic in- 
teractions, as well as the influence of curvature for small 
solutes. One may expect j s < 7, and lowering 7 S will lead 
to less hydrophobic attraction. A possible approach to in- 
clude the electrostatic field effects due to a net charge on 
the solutes is to absorb these effects into the solute-liquid 
surface tension. A naive procedure is to approximate 7 S 
by the solvation free energy per unit area of a charged so- 
lute. In section ITlI Bl it was shown that the Born expres- 
sion (4) fits the MD data for uniformly charged solutes 
quite well, once the Born radius Rb has been adjusted. 
Thus one may write for large solutes R > lnm (neglect- 



ing 1/e compared to 1): 



7s = 7 ' 



q z e z 



32ir 2 e R 3 R 



(14) 



which indeed lowers the hydrophobic attraction when 
charge is added to the solutes. No hydrophobic attraction 
occurs when 7 S = 0, so that the corresponding critical 
charge q c satisfies: 



32ir 2 e 2 e jR B 



(15) 



For Rb ~ lnm, one obtains \q c \ « 3, which in view of 
the sensitivity to Rb, yields at least the right order of 
magnitude, since the MD data discussed later suggest 
|g c | «8. 



C. Forces and potentials from simulation 

In the MD simulations, the mean force between two 
solutes was calculated by placing them at fixed positions 
Ri and R 2 along the body diagonal of the simulation 
cell, and averaging over water configurations generated 
during the runs which extended typically over 1-3 ns. 
The averaging was performed as long the statistical error 
was larger than A/3FA = 0.5, approximately twice the 
symbol size in the figures showing the forces. Separate 
MD simulations have to be carried out for each center-to- 
center distance Ri 2 = Ri — R 2 , i.e. for a series of surface- 
to-surface distances s = R\ 2 — 2R$. The force acting on 
solute 1 is estimated from the statistical average of the 
gradient of the total interaction energy V so \ in Eq. 

Fi(Ri2) = (-W so1 (R 12 ))b 12 : (16) 

where the constrained statistical average is taken over 
solvent configurations, when the two solutes are held 
fixed at a separation R\ 2 . Note that while the solute 
translational degrees are frozen, they rotate freely un- 
der the action of the torques exerted by the solvent and 
the other solute. In other words, the calculated effec- 
tive forces are orientationally averaged. By symmetry, 
^2(^12) = —Fi{Ri 2 ). The magnitude of the effective 
force is obtained by projecting onto the vector R\ 2 : 



R 
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(17) 



The resulting effective solute-solute potential follows 
from: 



w(R 12 ) = / F(R)dR. 



(18) 



In simulations where counterions are present, the sum 
of all pair interactions between the latter and the solute 
must be added to V^ i m Eq. (jl6|) . i.e. the mean force 
acting on the solute is the statistical average of the sum 
of the instantaneous forces exerted by all water molecules 
and ions on the solute, in the presence of a fixed second 
solute. 



V. MOLECULAR DYNAMICS RESULTS FOR 
PROFILES AND FORCES 

A. Neutral solutes 

We first consider the case of uncharged solutes. The 
water density profiles are illustrated in Figs. H3T aWe) 
for the case of solutes of radius R — llA and differ- 
ent surface-to-surface distances along the z-axis joining 
the centers. We plot density contours coded by vari- 
able shades of gray. The profiles are calculated using a 
cylindrical average around the symmetry axis (center-to- 
center line). The profiles show a considerable depletion 
(dark region) of the solvent between the spheres, reminis- 
cent of the observations of Wallquist and Berne for flatter 
solutes^ As the surface-to-surface distance s is increased 
for fixed radius R, the water molecules penetrate into 
the region between opposite solute surfaces, as signaled 
by a decreasing radius of the dark region between the 
spheres. Eventually at a distance between s f=s 5A and 
6 A (Figs. Etc) and (d)) the region between the solutes 
fills with liquid. When s > 7A, the solvent layers around 
an isolated solute are hardly disturbed by the presence 
of the other solute. 

Examples for the mean force for several radii 6A < R < 
13 A are shown in Fig. [21 The largest radii are of the or- 
der of the size of small globular proteins or of oil-in-water 
micelles. The average force obviously goes to zero at large 
distances s and for symmetry reasons, it is directed along 
the center-to-center axis. As expected from a depletion 
mechanism, the force is attractive and its contact values 
and range increase with R. We observe for all radii that 
the range of the force is similar to the distance at which 
the drying between the solutes vanishes. The potentials 
of mean force w(s) may be calculated for each R accord- 
ing to Eq. I|18|) . The resulting potentials are shown in 
the inset to Fig. They closely resemble results ob- 
tained for polymer-induced depletion potentials between 
spherical colloids, albeit on different length and energy 
scales^ The force at contact, F(s = 0), and the range of 
the force s c scales roughly with R in agreement with the 
macroscopic model prediction (14). From the MD data 
we extract a rough scaling F(s = 0) » 3.97-R and range 
s c w OAR, reasonably close to the theoretical estimates 
from Sec. IV. B. Overall there is a striking qualitative and 
even semi-quantitative agreement between the MD forces 
in Fig. 14 and the predictions of tbe macroscopic model 
in Fig. 12. 

B. Uniformly charged solutes 

We now turn to charged solutes. Consider first op- 
positely charged solutes. The water density profiles are 
illustrated in Figs. ll5f aWd1 for the case of solutes of ra- 
dius R — llA , a surface-to-surface distance s — 4A , and 
various charges q. The upper part of each frame shows 
a density contour plot coded by variable shades of gray. 
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FIG. 13: Contour density profiles of water around two 
neutral So solutes with radius R = llA for surface-to- 
surface distances (a) s — 2A, (b) s — 4A, (c) s = 5A, 
(d) s = 6A, and (e) s = 7A. The dark and light areas 
correspond to low and high densities of water molecules. 
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FIG. 14: Simulation results (symbols) of the mean force, 
F* = /3FA, between two neutral S solutes in SPC/E 
water. The lines are guides to the eye. Results are plot- 
ted for solute radii R = 6 A (squares), R — 9 A (di- 
amonds), R — 11 A (triangles pointing up), and R = 
13A (t riangles pointing left). The inset shows the inte- 
grated force (potential of mean force) in obvious order. 





The lower part shows density profiles along the center-to- 
center axis z, averaged over a coaxial cylindrical volume 
of radius 5A. In Frame (a) we plot the density distri- 
bution for the neutral case Q — for comparison with 
the charged systems. Frames (b)-(d) in Fig. ^| show 
water density profiles in the vicinity of two spheres car- 
rying opposite electric charges ±ge at their center (oppo- 
site charges ensure overall charge neutrality without any 
need for counterions). As q increases from zero (frame 
(a)), water is seen to penetrate between the two solutes, 
the central peak around z = in the density profiles 
increases rapidly and its amplitude reaches roughly the 
bulk density of water when q = 10. Note that this central 
peak is asymmetrically split, indicating the presence of 
two hydration layers which differ somewhat depending on 
their association with the anionic or cationic solute. This 
difference is also evident in the contact values of the out- 
side surfaces of the solutes, and is a consequence of the 
different arrangements of the water dipoles around the 
solutes induced by the local electric fields. The asym- 
metry of the profiles can be rationalized by inspecting 
the water structure around isolated solutes, as shown in 
Fig-EJa) and (b) and Fig. Urn Sec. IIII Al The hydration 
shell is more sharply defined around the cationic than 
around the anionic solute. The water dipoles tend to 
point radially away from the cation, while the opposite 
configuration is more favorable around anions. 

The resulting mean forces between solutes are plot- 
ted for q = 0, 2, 5 and 10, as functions of the surface-to- 
surface distance s in Fig. 1161 together with corresponding 
potentials of mean force. The mean force includes the di- 
rect Coulomb interaction between the two solutes (with 
proper account for the periodic images), which is in fact 



FIG. 15: Density profiles of the water molecules around 
(a) two neutral So solutes of radius R = llA and two op- 
positely charged S± solutes of radius R = 11 A carrying a 
charge ±qe with (b) \q\ = 2, (c) \q\ = 5, and (d) \q\ = 10. 
The surface-to-surface distance in all cases is s = 4A. 
In the contour plots dark regions indicate low density 
regions while high densities are plotted bright. The pan- 
els below the contour plots show the water density p 
scaled with water bulk density po in a cylinder of radius 
5 A, coaxial with the center-to-center line of the solutes. 



an order of magnitude larger than the total mean force. 
At large distances hydrophobic interactions become neg- 
ligible and the force should tend to — q 2 e 2 /(Aireoer 2 ), 
where r = 2R + s and e is the dielectric permittivity 
of bulk water; the corresponding curve for q = ±10 is 
also shown in Fig. 

The most striking result illustrated in Fig. is the 
near independence of the force at contact, s = 0, with 
respect to solute charge. From the density profiles in 
Fig. E| the hydrophobic attraction is expected to be re- 
duced but this reduction is almost exactly compensated 
by the Coulomb attraction between solutes in the pres- 
ence of the solvent. As q increases, the initial slope of 
the effective force increases. The potential of mean force 
(shown in the inset to Fig. I16|) exhibits a contact value 
which increases with q, indicating that the reduction of 
hydrophobic attractive free energy clearly outweighs the 
increase in bare Coulomb attraction between the latter. 
Simulations calculating the forces at and near contact 
for q = 7 and q = 15, not shown in Fig. 1161 confirm 
this trend. Note that the potential of mean force for 
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FIG. 16: Mean force, F* = (3Fk, as a function of 
surface-to-surface distance s for neutral So solutes (cir- 
cles) and oppositely charged S± solutes of R = llA with 
different central charges ±qe: q = 2 (squares), q — 
5 (diamonds), q = 10 (triangles up). The dashed 
line represents the electrostatic force between 2 peri- 
odically repeated solutes with opposite charges q = 
± 10 in a continuous solvent with permittivity e = 
80. The inset shows the resulting potentials of mean 
force; the contact values w(s = 0) increase with q. 




they are sharply peaked at a distance s « lA of the O 
atoms from the solute surface. This would lead to a com- 
plete shared hydration layer, and consequently to a kink 
in the force versus distance curve, between 2 flat solutes 
separated by s — 2A . This critical separation is shifted 
to shorter distances due to the curvature of spherical so- 
lutes. 

In view of this delicate balance between various in- 
teractions, we have also examined the case of equally 
charged solutes. In this case monovalent counterions 
(Na + or CI - ) were included to ensure overall charge neu- 
trality. The situation is summarized in Fig. ll7l for solutes 
of radius R = 11 A and charge q — ±4 and q = ±8. 
Charges of 4 and 8 were chosen to allow a direct compar- 
ison with the results for the models T± and C± in the 
next section. The water density profiles are symmetric 
with respect to z = for equally charged solutes, but 
differ substantially when going from a pair of anions to a 
pair of cations, as discussed in Sec. lIII Al This difference 
is reflected in the effective forces and potentials shown 
in Fig. El The interaction between the anionic solutes 
is always more attractive. For q = ±4 hydrophobic at- 
traction overcomes the repulsion between like charges, 
while for q ± 8 the electrostatic contribution dominates 
and the force is mainly repulsive, apart from a small at- 
tractive kink at s — 2 — 3A. The contact value of the 
repulsive forces is an order of magnitude higher than in a 
continuous solvent, originating obviously from the lack of 
water between the solutes and, hence a reduced dielectric 
screening. 

The reduction of the hydrophobic attraction between 
initially uncharged solutes, upon increasing the solute 
charge, may be qualitatively understood by the orient- 
ing action of the strong electric field between charged 
solutes on the water molecules, which disrupts the local 
hydrogen-bond structure and moves water locally away 
from conditions of liquid- vapor coexistence, so that "dry- 
ing" no longer occurs. 



C. Discrete charge patterns 



FIG. 17: Mean force, F* = /3FA, as a function of 
surface-to-surface distance s for equally charged S± so- 
lutes of radius R = 11 A and q = —4 (diamonds), 
q = 4 (circles), q = — 8 (squares), and q = 8 (triangles). 
The inset shows the according potential of mean force. 



q = 10 shows more long range attraction compared to 
the smaller q data due to the increased electrostatic at- 
traction. The eye-catching kink in the force for q =10, 
at a distance s ~ lA is reproducible, and is proba- 
bly related to the pronounced shell structure of water 
molecules around highly charged solutes, discussed in 
sec IIII Al While for neutral (and weakly charged) so- 
lutes, the O and H density profiles show little structure, 



The water density distribution around two tetrahe- 
dral solutes To (i.e carrying no net charge) is plotted 
in Fig. 18(a)-(e) for increasing values of the surface-to- 
surface distance s. As in the case of neutral solutes So, 
water depletion between the solutes is observable up to 
a solute distance of s ~ 5A. The bright regions near the 
surfaces indicate high density water and stem from the 
solvation shells in the immediate vicinity of the surface 
charges. We note again that the density profiles were 
calculated by averaging the water density cylindrically 
around the symmetry axis and the solutes were free to 
rotate in the simulations. The difference in brightness 
at the solute surfaces shows that certain orientational 
configurations are favored over others. If the tetrahedra 
were rotating freely (without any mutual interactions), 
the brightness would be the same everywhere on the so- 
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lute contour, as in the case of homogeneously charged 
solutes. It seems that the systems chooses those config- 
urations where the hydrophobic parts of the solutes (i.e. 
areas between the hydrophilic surface charges) face each 
other. We have also performed simulations of the overall 
neutral tetrahedral solute replacing the SPC/E water by 
a continuous solvent with permittivity e = 80. In the lat- 
ter case and for short distances, s < 5k, configurations 
are favored in which opposite surface charges associated 
with the two solutes face each other, thus strongly lower- 
ing the electrostatic energy of the system. With explicit 
water this is apparently no longer the case, despite the 
expected reduction in dielectric screening close to the sur- 
faces; the system tries to deplete water from the solute 
surfaces between the solutes. For distances larger than 
s ~ 5A , where no visible drying occurs, the positions 
of the bright regions are more smeared out on average, 
indicating a higher rotational freedom of the solutes. A 
simple configurational order parameter will be defined 
and discussed later. The water profiles for the tetrahe- 
dral solutes T + show the same behavior, in particular a 
visible drying for s < 5 A. We have not performed simula- 
tions of a pair of the negative T_ tetrahedra but expect 
similar behavior. 

The water density profiles around a solute Co carrying 
an overall neutral cubic charge distribution are plotted 
in Fig. ll9f aV(c1. No water depletion is visible at any dis- 
tance. In (a) the black region between the solutes comes 
from the fact that the solute surfaces touch. The posi- 
tions of the high density regions of water corresponding 
to the first solvation shells of the surface charges are on 
average distributed homogeneously over the sphere sur- 
face, pointing to a high orientational freedom of the so- 
lutes. The density of water close to the solute surface 
(bright ring in Figs. ll9f aWc11 is on average higher com- 
pared to the tetrahedra due to the larger surface charged 
density. The water density profiles for the positive cubic 
solute C+ show a different behavior, resembling the re- 
sults for the tetrahedral solute, as shown in Fig. [2UJ For 
distances s < 4A no bright region is found between the 
solutes, and hence water depletion is observed. Similar 
to the tetrahedral case the solutes stay mainly in orien- 
tational configurations in which the hydrophobic patches 
face each other. 

The effective force between the model solutes is plotted 
in Fig. |^] We show the force between pairs of neutral 
and overall positive tetrahedra, To and T + as well as 
between pairs of neutral and overall positive cubes, Co 
and C+. Simulations with overall charged solutes were 
carried out with explicit counterions. We have not per- 
formed simulations of a pair of overall negative tetrahe- 
dra and cubes. The force between pairs of overall neu- 
tral and charged tetrahedra is only slightly less attrac- 
tive than between pairs of spherically symmetric neu- 
tral solutes (compare to Fig. I14|l . This is surprising as 
one expects a stronger influence of the high electric fields 
generated by the surface charges. We learned from the 
investigation of homogeneously charged solutes that an 




FIG. 18: Contour density profiles of water around two 
To solutes for surface-to-surface distances (a) s = OA, 
(b) s = 2k, (c) s = 4A, (d) s = 6A, and (e) s = 9k. 
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FIG. 19: Contour density profiles of water 
around two Co solutes with cubic charge dis- 
tribution for surface-to-surface distances (a) 
s = OA , (b) s = 2k , and (c) s = Ak . 




FIG. 20: Contour density profiles of water around 
two C + solutes for surface-to-surface distances (a) 
s = OA , (b) s = 2k , and (c) s = 4A . 



electric field can considerably lower the hydrophobic at- 
traction. Apparently, a more anisotropic electric field 
distribution again favors hydrophobic attraction. Com- 
pare the force between pairs of positively charged tetra- 
hedra and pairs of homogeneously charged solutes with 
the same overall charge q = +4 (Fig. I17|) : the attrac- 
tion between the homogeneously charged solutes is less 




FIG. 21: Mean force between T+ (diamonds), 
T (triangles), C + (squares), and Co (circles) so- 
lutes. Also plotted is the force between two pe- 
riodically repeated C + solutes in a continuous sol- 
vent with e = 80 (solid line without symbols). 



than half that between T+ solutes. As already discussed 
above, the strong attraction between two tetrahedra is 
accompanied by water depletion between them, as in the 
case of homogeneous neutral solutes. The water density 
profile around an isolated tetrahedral solute T±, shown 
in Sec. IIII Al already illustrated strong depletion of water 
from the regions between the first solvation shells of the 
discrete surface charges. A possible explanation of the 
strong attraction between two tetrahedra is that this de- 
pletion is amplified when two hydrophobic patches of the 
solutes come close and face each other, and thus lowering 
the free energy. 

The effective force between two solutes with overall 
neutral cubic charge distribution shows qualitative dif- 
ferences compared to the tetrahedra. Cubes with zero 
overall charge still attract each other, but the interaction 
range is decreased. Analysis of the configurations shows 
that for close cubic solutes (s ~ 1 — 3A) one positive 
and one negative charge belonging to different solutes are 
on average very close, interacting with reduced dielectric 
screening than in bulk water due to their mutual prox- 
imity. The attraction observed is therefore mainly due 
to the electrostatic contribution and not hydrophobic at- 
traction as in the case of neutral tetrahedra. For the 
C+ solute the situation is again different. All charges re- 
pel, allowing the hydrophobic patches to face each other 
and water depletion is induced, as seen in the water pro- 
files of Fig. |20l Although equally charged, the cubic so- 
lutes still attract each other in striking contrast to the 
homogeneously charged solutes with Q = 8 in explicit 
water (Fig. I17|) and different to the case where water is 
replaced by a continuous solvent with e = 80, also plotted 
in Fig. ED 

In the following we investigate how the explicitly re- 
solved water affects the average orientational configura- 
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2R +s 




FIG. 22: Sketch of two T solutes with radius Rq 
and a surface-to-surface distance s. The dashed lines 
through the solute centers delimit a slab of width 
2R + s. The numbers of solute charges (small spheres) 
of each solute N\ and N2 inside the slab define the 
order parameter N — N± ■ N2, described in section 
IV CI N = 6 = 3 ■ 2 in the configuration shown. 




FIG. 23: Configurational probabilities P(N) as defined 
in Sec. IV CI for two non-interacting tetrahedral solutes. 



tions of two close tetrahedral solutes, when their centers 
are held at fixed positions. In a continuous solvent the 
probability of observing a certain configuration is purely 
determined by the electrostatic interactions between the 
surface charges. Obviously, for close distances of the so- 
lutes, structural effects of explicit water are expected to 
be very significant. A simple orientational order param- 
eter, which coarsely probes different orientational config- 
urations of a pair of tetrahedral solutes can be defined 
as follows: we count the number of surface charges in 



the slab delimited in width by the centers of the two 
solutes, as sketched in Fig. [221 Let N\ and N2 be the 
numbers of charges in the slab belonging to the first and 
second solutes. The values JVi, i=l,2 for a tetrahedron 
with four charged vertices are obviously Ni = 1,2 or 
3 (it is not possible to have four charges on one half 
sphere of the tetrahedral solute). The order parameter 
is now defined as the product N = N1N2 and can take 
values 1,2,3,4,6,9, which characterizes 6 different mutual 
orientational configurations. For N = 1, for instance, 
one charge of each solute is located in the slab, and the 
charges are both necessarily close to the symmetry axis; 
on the other hand, when N = 9, 3 charges of each so- 
lute are within the slab and the bare triangular surfaces 
between the three charges are mainly facing each other. 
In Fig. 122 we sketch two tetrahedra in a configuration 
N = 3 ■ 2 = 6. The probability distribution P(N) for 
two freely (without any interactions) rotating tetrahedra 
is plotted in Fig.EBl iV = 2 = 2-l = l-2, iV = 4 = 2-2, 
and iV = 6 = 2- 3 = 3- 2 are the most likely configura- 
tions, with probabilities P(2) ~ 0.21, P(4) =~ 0.4, and 
P(6) = 0.28. 

In Fig. !24f aWd1 we plot the probability distribution 
P(N) for interacting tetrahedra in a continuous solvent 
with permittivity e — 80. In Fig. 124( a) and (b) we show 
the result for a pair of overall neutral To solutes at dis- 
tances s = 3A and s = 9A . For the close distance 
the free rotator distribution is dramatically changed and 
the N = 1 and N = 2 configurations are strongly fa- 
vored. This is due to negative and positive charges from 
different solutes attracting each other at close distance. 
For the larger distance the electrostatic interactions are 
weaker and P(N) strongly resembles the free rotator dis- 
tribution again. In Fig. I24f c) and (d) we show the same 
distribution function, now for a pair of overall positive 
solutes at distances s = 3A and s = 9A , resp. Here, at 
close distance the N = 1 and N = 2 configurations are 
suppressed due to the proximity of like charges, and the 
N > 3 configurations are enhanced, since they allow the 
charges of one solute to be at larger mean distance from 
the like charges of the second solute. For large distances, 
(d), we again recover the free rotator distribution. 

In Fig. l25f a)-(d) the continuous solvent is now replaced 
by explicit water molecules. For close solutes (s = 3A) 
the difference with the continuous solvent is large: the 
probabilities of the N = 1 and N — 2 configurations 
are reduced both for the neutral (a) and overall charged 
tetrahedra (c). The N — 6 and N = 9 configurations 
are greatly enhanced, indicating that on average the are 
hydrophobic surfaces face each other, as expected from 
the water profiles in Fig. ll8f bl.fcl. Remarkably, even for 
case (a) the explicit solvent system strongly favors water 
depletion rather than the proximity of two unlike charges, 
which would lower the electrostatic energy significantly. 
Increasing the distance to s = 9A, the probabilities of 
the N = 6 and N = 9 configurations are lowered and the 
overall distributions, both for neutral and positive tetra- 
hedra are more similar to the free rotator distribution. 
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FIG. 24: Configurational probabilities P(N) as de- 
fined in Sec. IV CI for two tetrahedral solutes in a 
continuous solvent with e — 80. (a) and (b) are 
for T solutes at a surface-to-surface distance s = 
3A and s = Qk. (c) and (d) are for T + so- 
lutes at a distance s = 3A and s — 9A, resp. 
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FIG. 25: Same as in Fig. but now with ex- 

plicit SPC/E water instead of a continuous solvent. 



VI. CONCLUSION 

We have used a simple model of neutral and charged, 
nanometer-sized spherical solutes, embedded in explicit 
aqueous solvent, to investigate the influence of charge 
patterns on the solvation of a single solute, and on the ef- 
fective, solvent-induced interaction between two solutes. 
The charge patterns considered in this paper include uni- 
form charge distributions (equivalent to a single charge at 
the center of the spherical solute), as well as tetrahedral 
or cubic charge distributions, involving 4 or 8 discrete 
positive or negative charges situated at the solute sur- 
face, adding up to an overall positive, zero or negative 
charge Q (the T + ,T_,T and C+,C_,C models). 

Extensive constant pressure and constant temperature 



(NPT) Molecular Dynamics simulations were carried out 
under "normal" solvent conditions, i.e. close to liquid- 
vapor coexistence of water at room temperature. These 
simulations provide water density profiles around a sin- 
gle solute or a pair of solutes, which can be resolved into 
solute-oxygen and solute- hydrogen pair distribution func- 
tions, the distance resolved orientational order parame- 
ter P(r), the solvation free energy as a function of solute 
radius and charge, as well as the effective force and pair 
potential between two solutes, averaged over solvent con- 
figurations. The main results of this investigation may 
be summarized as follows: 

1. The density profiles of water around a single, neutral 
solute (So model), and their variation with solute radius 
R (cf. Fig. 3) exhibit the characteristic "destructuring" 
for radii R > 5A already reported by earlier studies4*§*£ 
The water molecules exhibit no significant orientational 
ordering around neutral nano-sized solutes. 

2. The hydrogen and oxygen density profiles change 
dramatically when the solute is uniformly charged (S± 
models). These profiles are sensitive to the anionic or 
cationic nature of the solute (for a given absolute charge 
|Q|), in addition the hydrogen profiles exhibit a splitting 
of the main peak in the case of anionic (S_) solutes (cf. 
Fig. 5). The orientational order parameter P(r) exhibits 
a significant structure, and a relatively slow decay with 
r, indicative of strong orientational ordering around the 
solutes S+ or S_, which is somewhat more pronounced 
around an anionic solute. The hydration asymmetry re- 
sults in preferential solvation of anionic solutes for a given 
radius and absolute charge \Q\, in agreement with earlier 
findings^*^ 

3. Moving from uniformly charged solutes to discrete 
(tetrahedral or cubic) charge patterns, the hydration of 
nano-sized solutes is found to exhibit a strong angular 
modulation associated with the hydrophilic "patches" 
around the discrete surface charges, and hydrophobic 
"patches" in between (cf. Fig. 6). The conflicting hy- 
dration patterns lead to a surprising depletion of water 
around T + or T_ solutes, compared to a neutral solute 
So- The solvation free energy is found to be about 20 % 
lower for solutes with discrete charge patterns compared 
to that of uniformly charged solutes with the same overall 
charge (cf. Fig. 9). 

4. The present simulations confirm the strong hy- 
drophobic attraction between two neutral spherical nano- 
sized solutes linked to solvent "drying", which was al- 
ready reported earlier for similar models. The MD results 
for solute radii R > 5A are nearly quantitatively repro- 
duced by a simple calculation based on purely macro- 
scopic considerations, and the force at solute-solute con- 
tact is found to scale roughly linearly with R. 

5. The effective attraction between neutral solutes is 
strongly reduced, or turns into a repulsion, when the 
nano-sized solutes carry equal, uniform charge distribu- 
tions. The total force has a repulsive electrostatic com- 
ponent, while examination of the water density profiles 
shows that the "drying" is mostly suppressed. The ef- 
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fective force is systematically less attractive (or more re- 
pulsive) between pairs of cationic solutes compared to 
anionic pairs (cf. Fig. 17). Turning to a pair of oppo- 
sitely (but uniformly) charged solutes, the MD simula- 
tions show that the range of the effective attraction de- 
creases when the absolute charge \Q\ increases, again in 
agreement with simple macroscopic considerations, but 
that the effective force at contact seems to be indepen- 
dent of Q, and equal to the hydrophobic force between 
neutral solutes; we have no explanation for this surprising 
observation. 

6. The situation for discrete solute charge patterns 
is, not surprisingly, more complex, due to the competi- 
tion between the resulting hydrophilic and hydrophobic 
"patches" on the solute surface. On average, some "dry- 
ing" of water is observed, and the resulting mean force 
between solutes carrying tetrahedral or cubic patterns 
is once more attractive, despite the electrostatic repul- 
sion ("like-charge attraction"). This effect is obviously 
incompatible with crude "implicit solvent" models. 

7. The complete break-down of "implicit solvent" 
models, whereby the latter is replaced by a dielectric 
continuum, is further illustrated by the highly coarse- 
grained representation of the configurational probability 
density of two solutes carrying discrete charge distribu- 
tions, introduced in Sec. V. The relative orientations of 
the surface charge patterns on the two solutes are com- 
pletely different for explicit and implicit solvent models, 
particularly at short surface-to-surface distances s (cf. 
Figs. 24 and 25). 

The key message of the present work is that explicit 
solvent models are unavoidable for a proper description of 
the effective interactions between nano-sized solutes like 
proteins, and that the latter are extremely sensitive to 
the precise location of any electric charges carried by the 
solutes. Contrarily to effective interactions on larger col- 
loidal scales, a generic coarse-graining strategy appears 
to be useless when solutes in the nanometer range are 
considered, and fully molecular models arc required for 
realistic simulations. 
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VIII. APPENDIX A: SIMULATION DETAILS 

The simulations were performed with the DLPOLY2 17 
package. The Berendsen barostat and thermostat were 
used to maintain the SPC/E water at a pressure of 1 
bar and a temperature T — 300K. For the simulations 
of the solutes with inhomogeneous charge distributions 
(T and C models) we used the rigid body algorithm with 



quaternions to properly account for the rotation of the 
anisotropic solutes. To this end we switched to an in- 
tegration routine using the Nose-Hoover barostat and 
thermostat, which turned out to be more stable in con- 
junction with quaternions. We carefully checked that 
both barostats give the same results by performing tests 
with bulk water, treated both with bond contraints and 
with the rigid body algorithm. Test runs using the Nose- 
Hoover barostat and thermostat for the S models also 
showed no difference. 

The simulation cell is a periodically repeated cube with 
a maximum boxlength of about L = 48A, containing up 
to iV w =3000 water molecules, depending on the solute 
size. For simulations of one isolated solute we required 
that the surface-to-surface distance to the nearest image 
solute was 2fjA, yielding a box size of L = 2R + 20A. 
For the calculation of the interaction force between two 
solutes, the latter are placed at fixed positions -Ri and 
R.2 on the body diagonal of the simulation cell. The cen- 
ter to center distance is then R12 = \Ri — i?2|- The 
corresponding box dimensions are chosen such that the 
surface-to-surface distance s = R12 — 2Rq to the near- 
est image solute is 20A. The box length can be calcu- 
lated as L = (4i? + s + 20A) / V3- Due to the constant 
pressure constraint the box length fluctuates slightly in 
the simulations. The long range electrostatic interac- 
tions were evaluated with smooth particle-mesh Ewald 
(SPME) summations 16 using 16 k vectors in each direc- 
tion and a convergence parameter of 3.2/r cut . A cutoff 
distance r cut = 9A was used for LJ-interactions and the 
real space SPME contributions. For the nanosized solutes 
a larger cut-off radius is obviously required. We optimize 
the computational speed by introducing a second cutoff 
for the solute- water and solute-ion interactions, chosen to 
be i?o + 4A, sufficient large for the shifted, short ranged 
repulsive interaction (Q. 



IX. APPENDIX B: FINITE CORRECTIONS 
FOR SOLVATION FREE ENERGIES FROM 
SIMULATION 

Accurate solavtion free energies for charged solutes 
can be obtained by Ewald summations in periodic cells 2 ^ 
when the self-interacion energy of the solute charges with 
its periodic images and the background charge is prop- 
erly included^ This correction is slightly modified when 
the solute size R is comparable to the box size The 
final expression for the electrostatic contribution to the 
solvation free energy for our solute models, including the 
finite size corrections, is: 

A M± = A/4 m 

e 2 e-1 /few v> 2 2^R 2 ^ 2 \ 

+ 8^~ — L4. + izrL9. 

\ a=l a—1 / 
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2 - 1 No Nc ( 1 \ 

a— 1 p^ta,l 

(19) 

where e is the macroscopic permittivity of water, and 
for a periodic array of cubic simulation cells, S,ew ~ 
—2.837297. In the case of a uniformly charged solute 
corresponding to a single charged site qe at the center, 
the result ifH^I reduces to2£ 



With the system sizes used in the present simulation the 
finite size corrections are very large and represent typi- 
cally twice the value of A^ s i m and stem mainly from the 
few-term in Eqs. (19) and (20). 



* e-mail address jd319@cam.ac.uk 

1 M. Rubinstein and R. H. Colby, Polymer Physics (Oxford 
University Press, 2003). 

2 D. Chandler (2004), to appear in Nature. 

3 C. N. Likos, Phys. Rep. 348, 267 (2001). 

4 F. H. Stilinger, J. Solution Chem. 2, 141 (1973). 

5 G. Hummer and S. Garde, Phys. Rev. Lett. 80, 4193 
(1998). 

6 K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 
103, 4570 (1999). 

7 D. M. Huang, P. L. Geissler, and D. Chandler, J. Phys. 
Chem. B 105, 6704 (2001). 

8 A. Wallquist and B. J. Berne, J. Phys. Chem. 99, 2893 
(1995). 

9 H. Shinto, M. Miyahara, and K. Higashitani, J. Coll. In- 
terface Sci. 209, 79 (1999). 

10 M. Kinoshita, S. Iba, K. Kuwamoto, and M. Harada, J. 
Chem. Phys. 105, 7177 (1996). 

11 Y. Qin and K. A. Fichthorn, J. Chem. Phys. 119, 9745 
(2003). 

12 E. Allahyarov, H. Lowen, A. Louis, and J. Hansen, Euro- 
phys. Lett. 57, 731 (2002). 

13 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. 
Phys. Chem. 91, 6269 (1987). 

14 J. Dzubiella and J.-P. Hansen, J. Chem. Phys. 119, 12049 
(2003). 

15 E. Spohr, Electrochim. Acta 44, 1697 (1999). 

16 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, 
H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 
(1995). 



17 W. Smith and T. R. Forester (1999), the DLPOLY.2 User 
Manual. 

18 D. Frenkel and B. Smit, Understanding Molecular Simu- 
lation: From Algorithms to Applications (Academic Press, 
1996). 

19 M. Born, Z. Phys. 1, 45 (1920). 

20 H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 
31, 369 (1959). 

21 W. M. Latimer, K. S. Pitzer, and C. M. Slansky, J. Chem. 
Phys. 7, 108 (1939). 

22 G. Hummer, L. Pratt, and A. E. Garcia, J. Phys. Chem. 
100, 1206 (1996). 

23 G. Hummer, L. Pratt, and A. E. Garcia, J. Chem. Phys. 
107, 9275 (1997). 

24 R. M. Lynden-Bell and J. C. Rasaiah, J. Chem. Phys. 107, 
1981 (1997). 

25 S. Rajamani, T. Ghosh, and S. Garde, J. Chem. Phys. 120, 
4457 (2004). 

26 P. G. Bolhuis and D. Chandler, J. Chem. Phys. 113, 8154 
(2000). 

27 A. A. Louis, P. G. Bolhuis, and J. P. Hansen, J. Chem. 
Phys. 117, 1893 (2002). 

28 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gun- 
steren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 
3684 (1984). 

29 J.-P. Hansen, Molecular Dynamics Simulation of Statistical 
Mechanical Systems (North Holland, Amsterdam, 1986), 
edited by G. Ciccotti and W.G. Hoover. 



